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Abstract 

We  consider  the  interrogation  by  means  of  a  pulsed  planar  electromagnetic  wave  of  a  dielectric 
slab  with  a  supraconductive  backing.  Previous  work  using  a  weak  formulation  with  finite  ele¬ 
ments  (FE)  demonstrated  the  ability  to  determine  material  parameters  and  the  slab  thickness  in 
the  inverse  problem.  In  this  work  we  report  on  results  using  Proper  Orthogonal  Decomposition 
(POD)  to  create  a  more  efficient  set  of  basis  functions  than  the  standard  FE  basis  functions. 
We  first  demonstrate  the  ability  of  the  reduced  basis  POD  formulation  to  capture  the  electro¬ 
magnetic  behavior  in  the  case  of  the  forward  problem.  We  then  apply  the  POD  formulation 
to  the  inverse  scattering  problem  with  unknown  parameters  and  show  that  the  POD  formula¬ 
tion  provides  a  considerable  reduction  in  computational  time  over  standard  FE  methods  with 
comparable  ability  to  recover  the  unknown  parameter  values. 
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1  INTRODUCTION 


Determination  of  material  properties  by  means  of  a  non-invasive  probes  such  as  low  energy  electromag¬ 
netic  pulses  is  desirable  in  a  broad  range  of  applications  in  biology  [1]  and  industry.  Previous  work  has 
shown  how  metallic  or  other  reflective  backings  could  be  exploited  to  obtain  information  about  mate¬ 
rial  properties  and  geometry  [10].  In  that  work,  a  finite  element  (FE)  based  variational  formulation 
was  employed  that  incorporated  Maxwell’s  full  equations,  the  antenna  source  current,  and  constitu¬ 
tive  polarization  models  (Debye  or  Lorentz).  Application  of  the  FE  variational  formulation  to  a  ID 
pulsed  input  scattering  problem  demonstrated  the  ability  to  capture  the  electromagnetic  behavior  in 
the  forward  problem.  Using  windowed  microwave  pulses  in  a  FE  variational  formulation,  the  authors  of 
[10]  were  able  to  effectively  utilize  the  two  reflected  signals  (from  the  air-dielectric  interface  and  from 
the  dielectric-backing  interface)  to  estimate  material  dielectric  properties  and  material  thickness  in  the 
inverse  problem. 

Solution  of  the  inverse  problem  to  identify  material  properties  can  be  time  consuming,  however,  since  it 
requires  repeated  calls  to  the  forward  finite  element  simulation,  whose  solution  time  grows  as  the  square 
of  the  number  of  nodes.  In  an  effort  to  reduce  computational  time  we  have  implemented  a  reduced  order 
formulation  of  the  variational  approach  of  Banks,  et  ai.[10].  In  place  of  the  finite  element  basis  functions 
our  variational  formulation  uses  basis  elements  obtained  through  application  of  the  method  of  proper 
orthogonal  decomposition  (POD)  to  standard  finite  element  computational  results. 

Our  efforts  reported  on  here  were  motivated  by  the  very  successful  use  of  POD-based  reduced  order 
methods  in  other  electromagnetic  interrogation  problems.  In  [9,  11]  the  authors  considered  eddy  current 
based  methods  for  interrogaton  of  dielectric  materials  for  anomalies  (flaws,  damages,  etc.).  Through 
both  computational  and  experimental  validation  efforts,  they  demonstrated  the  enormous  potential 
for  computational  savings  in  this  class  of  problems.  These  efforts,  however,  employed  eddy  currents 
originating  from  smooth  (uniform)  AC  signals  in  conductive  sheets.  This  permitted  the  reduction  of  the 
associated  time  domain  Maxwell’s  equation  system  to  a  phasor  form  involving  harmonic  systems  without 
transients.  A  question  of  great  interest  is  whether  the  reduced  order  ideas  of  [9,  11]  can  be  successfully 
employed  in  electromagnetic  interrogation  problems  where  the  interrrogating  signal  is  a  microwave  pulse 
and  the  resulting  dielectric  interface  reflections/transmissions  involve  transients  that  must  be  accurately 
computed  in  order  to  perform  the  material  identification  characterization  desired.  In  this  paper  we 
present  the  first  evidence  that  this  question  can  be  answered  in  the  affirmative. 

Proper  orthogonal  decomposition,  also  known  as  principal  component  analysis  [21]  and  Karhunen-Loeve 
expansion  [22,  26],  is  a  well  known  method  for  feature  extraction  in  statistical  and  pattern  recognition 
fields  [18].  The  POD  method  has  also  been  applied  in  a  wide  variety  of  other  fields  such  as  materials 
processing  [12,  23,  24,  29,  35],  characterization  of  human  faces  [33],  and  turbulent  coherent  flows  ([4,  14, 
15,  16,  19,  27,  34]  -  see  also  the  surveys  [13,  28]). 

The  POD  method  linearly  transforms  a  multivariate  data  set  into  an  optimal  set  of  uncorrelated  variables 
(POD  modes).  The  original  multivariate  data  can  be  written  as  linear  combinations  of  the  POD  modes. 
In  many  cases  the  POD  modes  more  efficiently  describe  the  variability  of  the  original  data  and  some 
dimensional  reduction  is  possible  by  retaining  only  the  most  important  modes. 

Recent  use  of  POD  for  reduction  of  order  in  distributed  parameter  systems  includes  as  noted  above 
applications  to  parameter  estimation  or  inverse  problems  [9,  11]  as  well  as  applications  to  both  open 
loop  and  feedback  control  design  [2,  3,  7,  8,  12,  23,  24,  25,  29,  30,  35].  Computational  evidence  from  a 
number  of  fluid  and  electromagnetic  applications  [12,  9,  11,  24,  23]  indicates  that  the  important  features 
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of  these  systems  are  essentially  of  low  finite  dimension  and  amenable  to  approximation  by  a  finite  span 
of  appropriately  chosen  basis  elements. 

We  apply  our  POD  reduced  order  model  to  a  ID  scattering  problem  in  which  a  windowed  microwave 
pulse  is  incident  on  a  dielectric  slab  with  a  supraconductive  backing.  We  present  results  for  the  forward 
problem  which  show  that  the  reduced  order  model  is  capable  of  capturing  the  electromagnetic  behavior 
using  significantly  fewer  modes  than  the  FEM.  We  then  apply  the  reduced  order  model  to  the  inverse 
problem  and  show  that  the  reduced  order  model  is  capable  of  accurately  estimating  the  parameters  with 
a  significant  reduction  in  computational  time  over  the  standard  FEM.  We  believe  that  this  work  demon¬ 
strates  the  feasibility  of  utilizing  the  POD  reduced  order  variational  method  for  parameter  estimation 
in  transient  electromagnetic  systems  when  computational  time  is  an  important  consideration. 


2  Problem  Formulation  and  Approximations 


The  physical  problem.  Consider  the  problem  of  interrogating  an  infinite  (in  the  x  and  y  directions) 
slab  of  homogeneous  material  by  a  windowed  microwave  pulse  (Fig.  1).  The  interrogating  signal  is 
chosen  to  be  a  polarized  planar  electromagnetic  wave  normally  incident  on  an  infinite  slab  of  material 
contained  in  the  interval  [zi,Z2]  with  faces  parallel  to  the  xy  plane.  The  electric  field  is  polarized  with 
oscillations  in  the  xz  plane  only.  The  slab  will  be  denoted  by  D  and  the  region  exterior  to  the  slab  by 
Do- 


x 


Figure  1:  Geometry  of  the  ID  scattering  problem. 


The  electric  and  magnetic  fields  in  fi  and  Do  are  governed  by  the  macroscopic  Maxwell’s  equations  (see, 
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for  example,  [20]) 

V  x  E 

V  x  H 
V -D 
V-P 

The  vector- valued  functions  E  and  H  represent  the  strengths  of  the  electric  and  magnetic  fields,  respec¬ 
tively,  while  D  and  B  are  the  electric  and  magnetic  fluxes,  respectively.  The  conduction  and  source 
current  densities  are  represented  by  Jc  and  Js,  respectively.  The  electric  and  magnetic  polarizabilities 
are  represented  by  P  and  M,  respectively.  The  electric  and  magnetic  susceptibilities  are  represented  by 
eo  and  po,  respectively  and  the  density  of  free  electric  charge  is  represented  by  p. 

The  domain  fio  is  treated  as  empty  space.  As  such,  M  =  0,  P  =  0,  and  Jc  =  0.  In  addition,  it  is 
assumed  that  all  parameters  necessary  for  the  determination  of  the  fields  are  known  in  this  region.  The 
source  current  Js  is  generated  by  an  infinite  antenna  along  the  x  axis  and  hence  will  only  be  nonzero  at 
these  points  in  flo;  its  time-varying  component  will  generate  the  interrogating  electromagnetic  waves. 


dD 
dt 
P , 

0. 


dB 

W’  ^ 
+  J j 


D 

B 

J 


eo  E  +  P, 
PoH  +  po  M, 
Jc  +  Js, 


(2.1) 


For  the  purposes  of  this  study,  we  will  make  certain  assumptions  about  the  material  properties  of  the 
domain  Q,  neglecting  magnetic  effects  and  assuming  that  Ohm’s  law  governs  the  electrical  conductivity. 
Thus  we  have 


M(x)  =  0 
Jc(x)  =  aE(x) 


(2.2) 


Polarization  model.  For  simplicity,  we  treat  the  instantaneous  polarization  to  be  proportional  to  the 
electric  field  so  that  pn  =  eo\E,  where  x  is  a  dielectric  constant.  The  electric  flux  density  in  (2.1)  then 
can  be  written  as 

D  =  eo  E  +  eoX-E  +  P 

=  e0(l+x)|  +  P  (2.3) 

=  t0erE  +  P , 


where  er  =  1  +  x  >  1  is  the  relative  permeability.  The  parameter  er  is  a  spatially-dependent  parameter 
that  allows  for  instantaneous  effects  on  the  the  displacement  in  fl  due  to  the  electric  field  originating  in 
flo-  The  remainder  of  the  macroscopic  electric  polarizability  of  the  medium  fl  is  denoted  by  P  and  is 
represented  as  an  integral  representation  dependent  upon  a  dielectric  response  function  (DRF)  and  the 
past  history  of  the  electric  field  (for  a  discussion  of  this  model  and  its  relationship  to  other  models  in 
the  literature,  see  [10],  pg.  10) 


P(t,x)=  /  g(t  —  s,x)E(s,x)ds, 
Jo 


(2.4) 


where  it  is  assumed  that  P(0,x)  =  0. 


In  particular,  we  choose  the  Debye  model  [10,  17]  for  orientational  or  dipolar  polarization  in  O  described 
by 

tP  +  P  =  eo(es  — eoo)P  /o  c\ 

D  =  e^eoE  +  P, 

where  es  is  the  static  relative  permeability.  The  Debye  model  corresponds  to  g{t)  =  e_t/,reo (eg  —  )/r 

in  (2.4)  and  models  the  behavior  of  materials  whose  molecules  possess  permanent  dipole  moments. 
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The  magnitude  of  the  polarization  P  represents  the  degree  of  alignment  of  the  individual  moments. 
Comparing  (2.5)  and  (2.3),  one  finds  that 


(  Coo  in  0 
(  1  in  fio- 


(2.6) 


for  the  Debye  polarization  model. 


Mathematical  formulation.  Using  (2.1)  and  the  condition  of  a  plane-polarized  interrogating  input 
signal,  one  can  argue  (see  [10])  that  E(t,x)  =  and  H(t,x)  =  jH(t,z )  for  scalars  E  and  H  in 

fio-  As  a  result,  the  polarizabilty  P  and  flux  density  D  are  confined  to  the  xz  plane  (with  nonzero 
components  only  in  the  x  direction)  and  can  be  denoted  by  scalar  values  P  and  D.  Since  the  material 
0  is  homogeneous  in  the  xy  direction,  the  fields  can  be  represented  by  scalar  values  E,  H,  P ,  and  D  in 
ff  as  well.  Recalling  that  Js  also  lies  in  the  x  direction  we  find  that  Maxwell’s  equations  then  reduce  to 


„  d  H 

7^ 


~dt 


+  (7  E  +  Js 


(2.7) 


By  taking  appropriate  derivatives,  we  can  combine  the  equations  (2.7)  to  give 

Ho  tE  +  HoP  +  Ho&P  —  E  =  —fjojsj  (2-8) 

where  (2.3)  was  used  to  eliminate  D,  e  =  eo(l  +  (e-r  —  1)/q),  Is  denotes  an  indicator  function  for  set  5 
( i.e .,  Is  =  1  in  5  and  0  otherwise),  '  =  Jj,  and  '=  Jj. 


Boundary  conditions.  An  absorbing  boundary  condition  is  placed  at  2  =  0  to  prevent  the  reflection 
(back  into  the  region  of  interest)  of  waves  so  that 


1  dE  dE' 
c  dt  dz 


(2.9) 


where  c2  =  1/eoHo-  The  boundary  conditions  at  the  slab-supraconductive  backing  boundary  are  B -n  =  0, 
which  is  automatically  satisfied,  since  B  =  jB ,  and  the  condition  Exh  =  0,  which  gives  \Eyi  —  Exj\r  =  0, 
so  that  Ex  =  Ev  =  0  on  the  boundary.  This  is  equivalent  to  E(t,  1)  =  0  in  our  system.  In  addition,  the 
source  current  Js  is  specified  at  z  =  0  (as  we  have  noted  along  the  x-axis). 


Method  of  mappings.  While  the  value  of  z±  is  known,  the  value  of  £2  is  generally  unknown.  We  use 
a  piece-wise  linear  mapping  that  leaves  the  interval  (0, 21)  invariant  and  maps  ( Z\,Z2 )  to  (zi,l).  The 
effect  of  this  mapping  [5,  6,  10,  31]  is  to  transform  the  original  unknown  geometry  fi0  U  ^  t°  one  with 
a  known  geometry  S7  =  [0, 1].  A  new  coordinate  variable  in  Cl  =  [0, 1]  can  be  defined  as 


5  =  f(z)  =  < 

\ 

0  <  z  <  zi, 

(2.10) 

[  Zl  Zl )  Z2  ~  Zi  ! 

Zi  <  Z  <  Z2  • 

This  can  also  be  written  as 

m 

=  z  +  {(  -  l){z  -  Zi)I[zi 

,22] 

(2.11) 

where  (  =  ^0_zz\  and  f'(z)  =  1  +  ((  —  l)In(z).  The  chain  rule  is  used  to  convert  spatial  derivatives  to 
derivatives  in  terms  of  the  new  variable  z.  That  is, 


d_ 

dz 


dzd_ 
dz  dz 


(2.12) 
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Using  dz  =  f'(z)dz,  one  can  write  the  inner  products  of  the  weak  form  as 

(<i>,  U)  =  [  <j>{. z)tp{z)dz 

Jo 

rl  (j>{z)i>{z) 


-  jfr 


+  (C  -  l)/n 

(<t ')  =  [  <t>'  (z)ip'  (z)dz 

J  0„ 


dz 


=  [  (1  +  (C  -  l)-fn)^'(2)V!'(i)di 

Jo 


(2.13) 


(2.14) 


Forward  or  simulation  problem.  Consider  the  case  of  a  Debye  medium  fl  with  a  supraconductive 
backing  and  er(z)  =  defining  the  instantaneous  polarization  in  Q.  We  first  express  (2.8)  in  an 
alternate  form  by  integrating  it  against  a  “test”  function  <f>,  obtaining 

(a*o eE,  +  ^0P,  </>) 

-  (E",  4)  =  ■),<!>), 

where  the  polarization  P  is  of  the  form  (2.4)  and  the  mapping  to  Cl  =  [0, 1]  has  already  been  carried  out. 


Integration  by  parts  of  (2.15)  gives  the  weak  form 

(a»o eE,  +  (no aE,  <^j  +  (no P,  <j)^ 

+  +  =  U"1(" 

where  the  term  \E{t,  0)^(0)  is  part  of  the  weak  form  resulting  from  the  absorbing  boundary  condition 
(2.9). 


For  computational  purposes  the  time  variable  is  scaled  by  a  factor  of  c  =  1/  ^/eo  no  (i  =  ct)  and  the 
polarization  P  by  a  factor  of  1/eo  (P  —  P/e o)-  We  assume  that  the  electric  permittivity  and  magnetic 
permeability  of  the  medium  0  are  constant.  The  scaled  equation  becomes 

(erE,  <t>^  +  rfy  (aE,(p^  + 

+  {E',4>')  +E(t,  0)0(0)  =  -Vo  (j.(t,  ■),<!>),  "'i,; 

where  er(z)  =  1  +  In(z)(eoa  —  1)  is  the  relative  electric  permittivity  so  that  e  =  ereo  and  the  impedance 
of  free  space  is  defined  r]o  =  sfnojeo  ~  376.73  Ohms.  The  inner  products  <  •  >  are  the  weighted  inner 

products  of  (2.13). 


FEM  discretization.  The  interval  [0, 1]  is  uniformly  divided  at  the  points  z1^  =  ih,  where  h  =  1/N 
and  i  =  0,  ...,1V.  The  electric  and  polarization  fields  E  and  P  are  discretized  with  N  piece- wise 
linear  spline  functions  4>f{z)  such  that  </>n*y)  =  tv  for  i.j  =  0,  ...,1V.  The  essential  boundary 
condition  <pP( 1)  =  0  is  satisfied  by  omitting  <f>^  in  the  finite  dimensional  approximating  subspace 
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VN  =  span{^,  <$ , . . . ,  i The  fields  are  then  approximated  by 


E(t,  z)  as  EN(t,  z)  =  ^  ef  (f)</>f  (z)> 

£r° 

P{t,z)  «  PN{t,z)  =  ^  P? {t)<p? {z) | 

i=0 


(2.18) 


In  order  to  avoid  computational  difficulties,  we  make  one  further  requirement  that  the  material  bound¬ 
aries  of  the  slab  [zi,  1]  coincide  with  grid  points.  The  grid  point  coinciding  with  the  material  boundary 
zi  will  be  denoted  by  j  =  L.  By  design,  the  right  edge  of  the  material  has  been  mapped  to  z  =  1,  which 
corresponds  to  the  grid  point  j  =  N. 


Substitution  of  the  approximations  (2.18)  into  (2.17)  yields 


iV-1 

(t)  {er<j>i  ,<t>) 

N- 1 

+  E 

Voe? (f)  (c T<t> f , 

<P)  +  E 

i= 0 

i= 0 
N- 1 

i= 0 
iV-1 

(2.19) 

+  E 

’  +  E  ef  Wf  (0)0(0)  =  - 

■Vo 

(< Js{t ,  -),<p)  , 

i= 0 

i  0 

which  can  be  written  as 

(Mn  +  M 

n  (eoo  -  1))^  +  (r)c<  +  BN)eN  +  M»pN  +  KN 

eN 

=  VoJN 

(2.20) 

for 

e  =  (eo  ■ 

PN  pN 

c2  j  •  •  •  j  ^N—: 

L)  and  p  = 

(Po  1P2  i  ■  ■  ■  j  Pn—i  )  •  The  elements  for  the 

N 

x  N  matrices 

are  given 

by 

—  (Iq  <Pi :  (f>j ) 

_  f  In<j>i4>j 

Jo  1  +  (C  -  1  )la  ’ 

Mfj 

=  {(pij  (pj  )  = 

f1  dz 

Jo  1  +  (C  -  l)Ia  ’ 

(2.21) 

n 

=  <PM4>M 

Kn- 

n 

=  {Mj)  = 

/  (1  +  (C  -  1  )In)4>i4>jdz, 
Jo 

while  the  N 

'  x  1  finite  element  vector  JN  is  given  by 

JiN  = 

f1  Mi 

Jo  1  +  (C-1  )IndZ' 

(2.22) 

Note  that  the  integrals  are  in  terms  of  the  scaled  variables  of  (2.10). 


We  will  use  the  Debye  polarization  model  (2.5)  to  provide  the  constitutive  law  relating  the  polarization 
P  to  the  electric  field  E.  Applying  the  same  time  scaling  as  above  (P  =  P/eo,t  =  ct)  we  obtain 

P  P  A P  —  td\E  in  ff,  (2.23) 

where  td,  =  es  —  and  A  =  1/cr.  To  generalize  the  above  equation  to  the  entire  domain  we  can  multiply 
the  equation  by  the  indicator  function  Iq(z).  Then,  applying  a  Galerkin  approximation,  we  obtain 

M£pN  +  M»\pN  -  Mq  \edeN  =  0. 


(2.24) 
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The  combined  system  of  equations  is 

(MN  +  Mq  (eoo  -  l))e-JV  +  MqPn  +  (rioaMg  +  B)eN  +  KNeN 

pN  +  A pN  —  \tdeN 

Substituting  the  second  equation  and  its  derivative  into  the  first,  one  finds 


=  Vo  JN 
=  0. 


(Mn  +  Mq  (too  —  l))eN  +  (\edM»  +  r,0aM»+B)eN 

+  (-\2edMg +  KN)eN +  \2MnpN  =  voJN 
pN  +  A pN  —  A  ed.eN  =  0. 

This  can  be  written  as  a  first-order  system  with  the  composite  variable  x  =  (eN ,pN ,eN)  as 


or 


where 


r  eN  1 

PN 

+ 

eN 

M3N  x  +  K3Nx 

=  f3N 

0N 

0  N 

-IN 

0  N 

M? 

m? 

MP 

\ eN  1 

-  QN  - 

PN 

= 

0  N 

eN 

VoJN 

Mf  =  Mn  +  (eoo  -  l)M^ , 
M?  =  -A  2edMg  +  KN, 

=  A  2M», 


Mf  = 


M£(\ed  +Voo)  +  B ", 


and  Ilr  is  the  TV  x  N  identity  matrix  where  the  ones  have  been  replaced  with  zeros  in  rows 

L  -  1. 


(2.25) 

(2.26) 

(2.27) 

(2.28) 

(2.29) 
1  through 


The  form  of  the  source  current  is  chosen  to  be 


Js(t,z)  =  8(z)gs(t)I(0,tf){t)  =  -<5(z)sin(wt)/(0,t/)(t),  (2.30) 

where  w  is  a  specified  angular  frequency  of  the  input  signal  (and  the  carrier  frequency  of  the  resulting 
planar  wave)  and  6(z)  is  the  Dirac  distribution  with  infinite  mass  at  z  =  0.  The  signal  is  truncated  at 
finite  time  t ;  coinciding  with  gs(t)  =  0  by  the  indicator  function  ( t )  to  avoid  complications  arising 

from  discontinuities  in  the  input  signal.  As  a  result,  ojtj  =  rnr  for  some  positive  even  integer  n.  The 
duration  of  the  pulse  must  be  sufficiently  short  to  distinguish  reflections  from  the  front  and  back  of  the 
medium,  thus  requiring  that  tfC  <  2(z2  —  Z\).  The  current  matrix  then  becomes 


JN 


U)COs{u)t)I(pti  )(f) 

o^-1 


(2.31) 


Solution  method.  We  use  the  standard  Crank-Nicholson  scheme  to  find  an  approximate  solution  for 
the  differential  equation  system  (2.27).  Rewriting  (2.27)  as 

x  =  f{t,  x)  =  (M3N)~1  (F3N  -  K3Nx)  (2.32) 

and  choosing  a  step  size  k.  we  make  the  iterative  approximation 


%n+ 1  —  %n  T  kfn-\-0 


(2.33) 


8 


where  xn  »  x(tn)  =  x(nk),  fn+e  =  (1  —  0)f(tn,xn)  +6f(tn+  i,xn+i),  and  for  the  Crank-Nicholson 
scheme  6  =1/2.  This  system  can  be  solved  directly  for  xn+i 


xn-\-i  —  xn  +  kyni  (2.34) 

where 

(M3N  +  k9K3N)yn  =  F3+0  -  R3Nxn ,  (2.35) 

Fn+e  =  (1  —  ®)FnN  +  ®Fn+ 1>  and  xo  =  0.  Equation  (2.35)  can  be  solved  by  means  of  block-Gaussian 
elimination  to  a  block  upper-triangular  system,  LU  factorization,  and  back-substitution. 

Construction  of  POD  modes.  Simulation  using  the  finite-element  system  (2.27)  above  provides  a 
multivariate  data  set  consisting  of  2 K  vectors 

X  =  {E?  ...E%  ...P?  ...P%}  (2.36) 

representing  the  N  nodal  values  of  E(t)  and  Pit)  at  K  time  points  during  the  simulation.  This  original 
data  set  X  is  transformed  to  a  new  set  of  uncorrelated  variables  (POD  modes)  by 

W  «&}=*$,  (2-37) 

where  the  columns  of  $  =  {</>fK,  4>\K , . . . ,  are  the  eigenvectors  of  the  product  matrix  (X'X)(f>3K  = 
A i<j>iK ,  ranked,  in  descending  order,  with  respect  to  the  associated  eigenvalue  and  the  prime  superscript 
denotes  the  transpose  of  the  matrix.  The  POD  modes  W  are  orthogonal  wf  ■  w ^  =  A iSij,  and  the 
transformation  of  variables  preserves  the  data  variability 

IK  IK  IK 

^2(X'X)kk  =  ^2(W'W)kk  =y£\k.  (2.38) 

k= 1  k=  1  k= 1 

Expansion  of  the  original  data  X  in  terms  of  the  most  significant  POD  modes  minimizes  the  mean  square 
error  of  a  reduced  basis  representation[18] 


(XR)ij='£(WU($ykj,  (2.39) 

k= 1 

where  R  <  2 K.  The  most  significant  POD  modes  are  those  corresponding  to  the  largest  eigenvalues, 
since  the  ratio  of  an  eigenvalue  to  the  summation  of  eigenvalues,  Xk/  J2f=i  -\b  giyes  the  percentage  of 
the  mean  square  error  unaccounted  for  by  eliminating  the  corresponding  POD  mode  w £  in  the  reduced 
basis  representation  [18].  The  best  stopping  point  in  the  expansion  (2.39)  depends,  in  general,  upon  the 
application  and  various  algorithms  have  been  proposed[21]. 

POD  discretization.  A  second  discretization  formulation  produces  the  reduced  basis  model.  In  this 
case,  we  first  use  the  POD  modes  { w % }  to  obtain  the  POD  elements 

JV 

^k(z)  =  ^2(wk)iipi(z),  k  =  1,2,...,2K  (2.40) 

i= 1 

where  the  functions  ipi(z )  are  the  finite  element  linear  interpolation  functions.  The  electric  and  polar¬ 
ization  fields  are  approximated  as  linear  combinations  of  the  POD  basis  elements  corresponding  to  the 
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most  significant  POD  modes 


E(t,  z)  fa  ER(t,z) 
P(t,z)  fa  PR(t,z ) 


i= 1 
R 


«=1 


(2.41) 


where,  in  this  case,  R  <  N.  Application  of  this  approximation  in  (2.17)  (in  this  case  we  use  POD  test 
functions  <f» ,  =  H/j  i  =  1, 2, . . . ,  R)  yields  a  Galerkin-POD  system 


R 


^ef(f)  (erVR,yR) 


which  can  be  written  as 


+  ef(t)<CT«f,«fl)+^pf(t)<*f,*fl) 

i= 1  i=l 

+  £  ef  (t)  <*?' ’,  *"')  +  Y,  e?m?mR(0)  =  (Ut,  ■),  *")  , 

i=l  i= 1 

(2.42) 


(Mfl  +  (£oo  -  l))efl  +  {r)o(?Mg  +  BR)eR  +  MgpR  +  KReR  =  Vo  JR  (2.43) 

for  e  =  (ef ,  ef , . . . ,  eR ),  p  =  (pR,pR,  ■  ■  ■  ),  and  e  =  ef ,  ef , . . . ,  eR ).  The  elements  for  the  R  x  R 


matrices  are  given  by 

Mgj 

_  P  la§*§* 

Jo  l  +  (C-l)/n  ’ 

Mg 

II 

a;.» 

II 

r  1  d/fivirfj 

/  *  J  dz 

Jo  1  +  (C-1  )la  ’ 

(2.44) 

Bg 

=  «f(0)«f(0) 

Kjj  =  <4>f ,  )  =  Al  +  (C  -  l)/n)«f'*f  dz, 

Jo 

while  the  fix  1  finite  element  vector  is 


<2-45> 

and  the  integrals  are  in  terms  of  the  scaled  variables  defined  in  (2.10).  Choosing  the  source  current  as 
in  (2.30)  we  find 

JR  =  ucos{u:t)Ii0,tf)(t)^R  (0).  (2.46) 

Application  of  the  POD  approximations  to  the  Debye  constitutive  law  (2.5)  relating  the  polarization  P 
to  the  electric  field  E  gives 

Mg pR  +  Mg\pR  -  Mg XedeR  =  0.  (2.47) 

Since  the  POD  basis  functions  Uf#  are  global  spatially  distributed  functions,  the  matrix  Mg  is  not 
singular,  which  differs  from  the  finite  element  formulation. 


The  combined  system  of  equations  is  (2.43)  and  (2.47)  is 

(. MR  +  Mg  {too  -  l))eR  +  MgpR  +  {rioaMg  +  BR)eR  +  KReR  =  r}0JR 

pR  +  X pR  -  \edeR  =  0. 


(2.48) 


10 


Substituting  the  second  equation  and  its  derivative  into  the  first,  we  obtain 


(MR  +  MR(e  oo-l))^ 


+  (A  edMR  +  r)0aMR  +  BR)eR 
+  (-X2edMg  +  KR)eR  +  X2MgpR  =  r]0JR 
IpR  +  A IpR  -  A edIeR  =  0, 


(2.49) 


where  I  is  the  R  x  R  identity.  This  can  be  written  as  a  first-order  system  with  the  3 R  vector  variable 
x  =  ( eR,pR ,  eR)  as 

M3Rx  +  K3Rx  =  F3R  (2.50) 


where 


r  eR  1 

0R 

0R 

-IR  ' 

f 1 

0R 

pFL 

+ 

-A  edI 

XI 

0R 

PR 

= 

0R 

.  eR . 

L  MR 

mr 

% 

.  ZR . 

1 

o 

Sr 

_ i 

Mr  =  MR  +  (eao-l)MR, 
MR  =  -A  2edMR  +  KR, 

MR  =  A  2MR, 

MR  =  MgiXea+vo^+B11. 


(2.51) 


(2.52) 


The  solution  to  the  differential  equation  system  (2.50)  can  be  approximated  using  a  Crank-Nicholson 
scheme  in  the  same  manner  as  described  above  for  the  FEM  differential  equation  system  (2.27). 

Inverse  problem.  We  formulate  the  problem  described  in  detail  in  [10]  of  the  interrogation  of  a 
Debye  dielectric  medium  with  a  supraconductive  backing  by  a  plane-polarized  windowed  wave.  In  the 
inverse  problem  we  attempt  to  determine  parameter  values  of  the  dielectric  material  from  experimental 
measurements  at  z  =  0  of  the  reflected  electric  field  .  This  is  accomplished  by  minimizing  the  L 2 
difference  between  the  experimental  data  and  simulation  results 

s 

J(q)  =  E  I  °;  $)  -  E*  I2’  (2-53) 

i—1 

where  5  is  the  number  of  sample  data  points  Ei,  at  uniform  time  intervals  U,  Q  is  the  set  of  admissable 
parameters,  and  E(ti,0;q)  are  the  electic  field  values  arising  from  simulations  with  parameters  q.  The 
vector  q  typically  contains  dielectric  and/or  conductivity  parameters  characterizing  the  medium  and/or 
geometric  properties  (see  [10]). 

For  our  computational  testing  of  the  methods  and  algorithms,  synthetic  data  £)  is  produced  by  adding 
random  noise  to  the  results  of  FEM  simulations  with  a  known  set  of  parameters,  i.e., 

Ei  =  Ei(l  +  vr]i),  (2-54) 

where  Ei  are  the  values  sampled  from  the  solution  with  known  parameters,  ip  are  independent  normally 
distributed  random  variables  with  mean  zero  and  variance  one  [10].  The  amplitude  of  the  noise  is 
proportional  to  the  signal  level  Ei  and  the  coefficient  v.  By  choice  of  v,  we  can  control  the  relative  noise. 
For  example,  10%  relative  noise  is  achieved  (with  probability  0.9545)  by  choice  of  v  —  0.05  since  rp  lies 
in  [—2(7,  2a\  =  [—2,2]  (with  probability  0.9545). 

Minimization  of  J(q)  is  performed  using  a  Broyden-Fletcher-Goldfarg-Shanno  (BFGS)  variable  metric 
algorithm  ( dfpmin( ),  pg.  428  [32]).  The  gradient  of  J(q)  is  obtained  from  a  forward-difference  approx¬ 
imation  using  an  algorithm  adapted  from  ( fdjac( ),  pg.  388  [32]).  The  BFGS  algorithm  and  related 
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functions  are  modified  for  double-precision.  In  addition,  constraints  are  added  to  reflect  the  physical 
limitation  of  parameters,  ie..  no  negative  values.  The  convergence  requirement  for  zeroing  the  gradient 
is  1  x  10-4  for  all  the  minimization  problems  reported  here. 

The  values  E(ti,0;q)  appearing  in  (2.53)  are  obtained  from  simulations  using  either  FE  or  POD  basis 
elements.  Simulation  using  FE  basis  vectors  is  straight-forward,  but  slow,  since  the  basis  vectors  are 
independent  of  the  parameter  values.  Minimization  of  J{q)  using  simulations  with  POD  basis  elements  is 
faster,  but  more  complicated,  since  the  POD  basis  elements  are  generated  from  snapshots  obtained  from 
FE  simulations  with  specified  parameter  values.  For  the  minimization  problem  a  collection  of  snapshots 
from  simulations  covering  a  range  of  parameter  values  ( q ))  are  used  to  generate  the  POD  basis  elements. 

XnK  =  {^2if(9l)  ■  ■  -X-2. 'k(Qti)},  (2.55) 

where  each  snapshot  represent  N  nodal  values  and  there  are  K  snapshots  each  for  the  electric  and 
polarization  fields  per  simulation.  In  general,  the  simulation  results  change  gradually  with  changes  in 
parameter  values  and  the  POD  modes  are  still  able  to  efficiently  represent  the  range  of  data. 
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3  COMPUTATIONAL  RESULTS 


Forward  problems.  We  simulated  the  interrogation  of  a  dielectric  slab  with  a  forward  face  at  z\  =  0.33 
and  a  supraconductive  backing  at  z-i  —  0.89  by  a  plane-polarized  windowed  wave  using  the  techniques 
outlined  above.  The  mapped  region  0  was  discretized  with  N  =  450  and  L  =  150.  The  dielectric 
parameter  values  were  es  =  35,  =  5,  r  =  1  x  10-11  s,  a  =  1.0  x  10-2  Ohm-1.  The  signal  parameters 

were  chosen  to  be:  <j>  =  1.8  x  109  s-1,  u  =  2n  x  <j>  rad/sec,  and  tf  =  ^  =  3.33  x  10-9  s.  The  differential 
system  was  integrated  with  a  step  size  of  dt  =  1.0  x  10-4  ns  for  a  total  of  =  10  ns. 

We  experimented  with  different  snapshot  intervals  for  the  creation  of  the  POD  modes.  POD  modes 
created  using  a  sample  rate,  or  snapshot  interval,  of  0.005  ns  provided  an  efficient  representation  of  the 
data.  Figure  2  shows  the  percent  variability  captured  for  the  field  data  (E(t)  and  P(t))  as  a  function 
of  the  number  of  modes  at  the  0.005  ns  sampling  interval.  Some  representative  values  are  99.099%, 
99.962%,  and  99.999%  for  30,  50,  and  72  modes,  respectively. 


Figure  2:  Percent  of  field  (E(t)  and  P(t))  variability  represented  as  a  function  of  the  number  of  modes 
used.  The  sampling  interval  was  0.005  ns,  giving  2001  snapshots  for  each  field. 
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We  illustrate  the  efficacy  of  the  POD  reduced  order  representation  by  simulating  the  forward  problem 
using  the  first  72  POD  modes  (of  4002  possible  modes)  generated  from  snapshots  taken  at  the  0.005  ns 
sampling  rate. 

Figures  3-6  compare  electric  field  results  obtained  using  FEM  (N=450  and  L=150)  and  POD  (R=72) 
simulations.  In  Figure  3  the  electric  field  is  plotted  as  a  function  of  the  mapped  distance  z  at  0.7  ns.  At 
this  time  the  windowed  pulse  has  not  yet  reached  the  dielectric  material.  Figure  4  depicts  the  electric 
field  at  5.0  ns,  where  there  are  both  reflection  and  transmission  pulses  originating  from  the  initial  pulse 
at  z=0.  In  both  cases  there  is  good  agreement  between  the  FEM  and  POD  simulations. 


Figure  3:  The  electric  field  at  0.7  ns  as  a  function  of  2  (mapped)  for  the  FEM  simulation  (solid  line) 
and  POD  simulation  (dash-dot  line).  The  POD  reduced  order  model  used  70  modes. 
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Figure  4:  The  electric  field  at  5.0  ns  as  a  function  of  2  (mapped)  for  the  FEM  simulation  (solid  line) 
and  POD  simulation  (dash-dot  line).  The  POD  reduced  order  model  used  70  modes. 
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Figures  5  and  6  show  the  electric  field  at  7.0  and  10.0  ns.  At  these  times,  only  the  transmitted  portion 
of  the  initial  pulse  remains.  The  discrepancies  between  the  FEM  and  POD  simulations  are  more  evident 
due  to  an  increase  in  relative  error  because  of  the  very  small  magnitude  of  the  electric  field  outside  the 
dielectric  material.  The  Brillouin  precursors  are  evident  from  the  large  amplitudes  of  the  forward  and 
trailing  peaks  in  both  figures.  There  is  still  excellent  agreement  between  the  FEM  and  POD  simulations. 


Figure  5:  The  electric  field  at  7.0  ns  as  a  function  of  z  (mapped)  for  the  FEM  simulation  (solid  line) 
and  POD  simulation  (dash-dot  line).  The  POD  reduced  order  model  used  70  modes. 
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The  solution  time  with  the  POD  model  (5  min)  is  significantly  less  than  the  solution  time  required  for 
the  450  node  FE  model  (272  min).  Since  the  efficiency  of  the  modes  is  relatively  independent  of  the 
number  of  nodes  used,  the  time  savings  of  the  POD  reduced  order  model  over  the  FE  model  becomes 
more  significant  when  more  nodes  are  required. 

This  point  is  illustrated  in  Figure  7,  where  it  can  be  seen  that  the  FEM  solution  time  is  proportional  to 
the  number  of  nodes  squared.  Also  plotted  are  the  equivalent  POD  solution  times,  which  are  calculated 
in  the  same  manner  as  above.  That  is,  for  each  N- noded  FEM  simulation  (the  parameter  values  are  the 
same  as  above,  except  =  5.0  ns)  the  POD  method  is  applied  to  snapshots  at  0.005  ns  intervals. 
The  simulation  time  using  the  POD  formulation  with  M  modes  (chosen  to  capture  at  least  99.999%  of 
the  variability)  is  the  equivalent  POD  solution  time.  Figure  7  shows  that  the  equivalent  POD  solution 
times  are  relatively  independent  of  the  number  of  FEM  nodes. 


N2  xIO5 

Figure  7:  The  solution  time  as  a  function  of  the  number  of  nodes  squared  for  the  FEM  (o)  with  tfinai= 5.0 
ns.  Also  shown  are  the  equivalent  POD  solution  times  (x),  where  the  number  of  modes  captures  at  least 
99.999%  of  the  data  variability.  In  both  cases,  data  points  are  connected  by  straight  lines. 
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Inverse  problem:  Example  I.  We  turn  to  computational  findings  for  a  one  parameter  (es)  problem 
of  the  interrogation  of  a  dielectric  slab  with  a  supraconductive  backing  by  a  plane-polarized  windowed 
wave.  In  the  inverse  problem  we  attempt  to  determine  the  parameter  values  (es)  of  the  dielectric  material 
from  measurements  of  the  electric  field  at  z  =  0.  The  dielectric  slab  dimensions  and  parameters  are  the 
same  as  above  in  the  forward  problem,  except  for  N  =  150,  L  =  50,  and  t/  =  |  =  1.11  x  10-9. 

In  this  example,  we  investigate  the  inverse  problem  for  the  single  unknown  parameter  es  in  the  range 
31  <  es  <  39.  Three  FE  simulations  (es  =  31,33,  and  39)  each  with  K  =  501  snapshots  of  the  N  —  1 
nodal  values  for  the  electic  and  polarization  fields  are  used  to  obtain  the  POD  basis  elements  for  the 
minimization  problem 

X*Kl  =  {X^K1^  =  31)  X^-1(es  =  35)  X^~\es  =  39)}.  (3.56) 


The  elecric  field  is  measured  at  0.01  ns  intervals  from  0  <  t  <  5  ns.  Figure  8  shows  the  observed  values 
of  the  electric  field  at  z  =  0  as  a  function  of  time. 


Figure  8:  The  observed  electric  field  at  z=0  as  a  function  of  time.  Measurements  are  made  at  0.01  ns 
intervals  from  0  to  5  ns. 
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The  ability  of  the  POD  formulation  to  recover  the  true  values  in  the  inverse  problem  is  investigated 
as  a  function  of  the  number  of  modes  used  (Fig.  9).  The  following  conditions  were  used  to  obtain 
these  results:  0%  noise,  es(true)=33,  es  (initial)  =37,  and  J(q)  minimized  with  respect  to  observation 
measurements  from  times  0  <  t  <  3.5  ns.  As  expected,  the  relative  error  generally  decreases  as  the 
number  of  modes  increases.  Oscillations  may  be  due  to  the  influence  of  the  addition  of  particular  modes 
with  respect  to  the  true  parameter  value.  As  Fig.  9  indicates,  the  relative  error  decreases  slowly  after 
approximately  50  modes. 


Figure  9:  The  percent  relative  recovered  parameter  error  (100|es(true)  -es (recovered) )/eg (true) |)  as  a 
function  of  the  number  of  modes  (circles).  The  curve  drawn  through  the  data  is  a  cubic  smoothing 
spline  ( csaps  in  Matlab)  with  smoothing  parameter  equal  to  0.1 
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We  next  investigated  the  ability  of  the  POD  formulation  to  recover  es  with  a  noisy  signal  and  a  fixed 
number  of  modes  (60).  Figure  10  gives  an  example  of  the  observed  electric  field  without  noise  and  with 
5%  noise  added.  In  this  case  es  =  33. 


Example  of  observed  and  noisy  signal  at  5%  error 


Figure  10:  A  portion  of  the  observed  electric  field  at  z=0  with  5%  noise  (solid  line)  and  without  added 
noise  (dashed  line),  for  es(true)  =  33. 
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The  accuracy  of  the  POD  formulation  in  the  above  inverse  problem  with  one  unknown  parameter  (es) 
is  compared  to  the  FEM  formulation  in  Fig.  11,  where  the  relative  error  is  plotted  as  a  function  of  the 
percent  noise  in  the  signal.  For  basis  of  comparison,  the  same  seed  (-123456)  for  the  random  number 
generator  is  used  for  all  cases  with  non-zero  noise.  In  general  the  relative  error  increases  linearly  with 
respect  to  the  percent  noise  for  both  the  POD  and  FEM  formulations.  Except  at  the  zero  percent  noise 
level,  the  POD  method  is  slightly  more  accurate  than  the  FEM  formulations.  We  attribute  this  to 
serendipity  rather  than  any  remarkable  methodological  principle. 


Figure  11:  The  percent  relative  recovered  parameter  error  (100|es(true)  -es (recovered) )/eg (true) |)  as  a 
function  of  the  percent  noise  for  the  POD  (o,  60  modes)  and  FEM  (x)  formulations.  In  both  case  data 
points  are  connected  by  straight  lines. 


A  plot  of  the  ratio  of  FEM  solution  times  to  POD  solution  times  illustrates  the  time  savings  offered  by 
the  POD  method.  As  Fig.  12  indicates,  the  POD  method  is  most  efficient  at  the  0%  noise  level  where  it 
is  nine  times  faster  than  the  FEM  formulation.  As  the  noise  level  increases  the  POD  efficiency  decreases, 
eventually  reaching  a  5.4-fold  time  saving  at  the  15%  noise  level.  The  POD  times  are  not  changing  much 
as  the  noise  level  changes,  but  the  FEM  solution  times  decrease  somewhat. 
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Figure  12:  The  ratio  of  solution  times  (Tfem/^pod)  as  a  function  of  the  percent  noise.  The  data 
points  are  connected  by  straight  lines. 


Inverse  problem:  Example  II.  In  this  example,  we  investigate  the  inverse  problem  with  r  as  the  single 
unknown  parameter.  The  problem  conditions  are  the  same  as  in  the  previous  example,  except  that  now 
we  employ  fixed  es  =  35  and  r  is  varied.  Three  FE  simulations  (r  =  1  x  10-12, 1  x  10-11  and  1  x  10-10) 
each  with  K  =  501  snapshots  of  the  N  —  1  nodal  values  for  the  electic  and  polarization  fields  are  used 
to  obtain  the  POD  basis  elements  for  the  minimization  problem 

X*Kl  =  =  1  x  IQ”12)  X?~\t  =  1  x  IQ”11)  X?~\r  =  1  x  10-10)}.  (3.57) 


The  ability  of  the  POD  formulation  to  recover  the  true  values  in  the  inverse  problem  is  plotted  as  a 
function  of  the  number  of  modes  used  in  Fig.  13.  The  following  conditions  were  used  in  these  results: 
0%  noise,  r(true)=5  x  10— 11 ,  r(initial)=3.5  x  10— 11 ,  and  J(q)  minimized  with  respect  to  observation 
measurements  from  times  0  <  t  <  3.5  ns.  The  relative  error  generally  decreases  as  the  number  of 
modes  increases.  As  in  the  previous  example,  oscillations  may  be  due  to  the  influence  of  the  addition 
of  particular  modes  with  respect  to  the  true  parameter  value.  The  relative  error  decreases  slowly  after 
approximately  40  modes,  where  it  is  already  under  1%. 
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NUMBER  OF  MODES 

Figure  13:  The  percent  relative  recovered  parameter  error  (100|r(true)  -r(recovered))/r(true)|)  as  a 
function  of  the  number  of  modes  (circles).  The  curve  drawn  through  the  data  is  a  cubic  smoothing 
spline  ( csaps  in  Matlab)  with  smoothing  parameter  equal  to  0.1 


The  accuracy  of  the  POD  formulation  with  r  unknown  is  compared  to  the  FEM  formulation  in  Fig.  14, 
where  the  relative  error  is  plotted  as  a  function  of  the  percent  noise  in  the  signal.  For  basis  of  comparison, 
the  same  seed  (-123456)  for  the  random  number  generator  is  used  for  all  cases  with  non-zero  noise.  As 
expected,  the  relative  error  increases  with  respect  to  the  percent  noise  for  both  the  POD  and  FEM 
formulations.  The  FEM  method  is  more  accurate  than  the  POD  at  low  noise  levels.  However,  above  the 
5%  noise  level  the  POD  method  is  increasingly  slightly  more  accurate  than  the  FE  method. 
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Figure  14:  The  percent  relative  recovered  parameter  error  (100|r(true)  -r(recovered))/r(true)|)  as  a 
function  of  the  percent  noise  for  the  POD  (o,  60  modes)  and  FEM  (x)  formulations.  In  both  case  data 
points  are  connected  by  straight  lines. 


A  plot  of  the  ratio  of  FEM  solution  times  to  POD  solution  times  illustrates  the  time  savings  offered  by 
the  POD  method.  As  Fig.  15  indicates,  the  POD  method  is  4.5  to  8  times  more  efficient  than  the  FEM 
formulation.  Variations  in  the  efficiency  are  largely  due  to  variations  in  the  solution  times  for  the  FEM. 
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Figure  15:  The  ratio  of  solution  times  (Tfem/^pod)  as  a  function  of  the  percent  noise.  The  data 
points  are  connected  by  straight  lines. 


Inverse  problem:  Example  III.  In  this  example  both  parameters  es  and  r  are  simultaneously  sought 
in  a  two  parameter  inverse  problem.  Except  for  es  and  t,  simulation  conditions  were  the  same  as  in  in  the 
one  parameter  inverse  problems  above.  The  elecric  field  is  measured  at  0.01  ns  intervals  from  0  <  t  <  5 
ns.  Nine  FE  simulations  (combinations  of  es  =  31,33,39  and  r  =  1  x  10-12, 1  x  10— 11 , 1  x  10-10)  each 
with  K  =  501  snapshots  of  the  N  =  149  nodal  values  for  the  electic  and  polarization  fields  are  used  to 
obtain  the  POD  basis  elements  for  the  minimization  problem. 

Figures  16  and  17  summarize  the  two  parameter  inverse  results  under  the  following  conditions  es(true)=33, 
es (initial)  =37,  r(true)=5  x  10— 11 ,  r(initial)=3.5  x  10— 11 ,  and  J(q)  minimized  with  respect  to  observation 
measurements  from  times  0  <  t  <  3.5  ns. 

Figure  16  plots  the  relative  error  for  the  recovered  es  values  as  a  function  of  the  noise  level  (0-20%)  for 
the  POD  and  FEM  methods.  From  Fig.  16  it  can  be  seen  that  the  FEM  results  are  generally  more 
accurate  than  the  POD  results,  but  that  the  relative  error  for  both  is  under  0.2%.  The  relationship 
between  the  relative  error  and  the  noise  level  appears  to  be  nonlinear,  unlike  the  one  parameter  recovery 
of  es. 
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Figure  16:  The  percent  relative  recovered  parameter  error  (100|es(true)  -es (recovered) )/es (true) |)  as  a 
function  of  the  percent  noise  for  the  POD  (o,  50  modes)  and  FEM  (x)  formulations.  Data  points  are 
connected  by  straight  lines. 


Figure  17  plots  the  relative  error  for  the  recovered  r  values  as  a  function  of  the  noise  level  (0-20%)  for 
the  POD  and  FEM  methods.  The  relative  error  appears  to  be  roughly  linear  with  respect  to  the  noise 
level  for  both  POD  and  FEM  formulations.  The  relative  error  for  both  formulations  do  not  differ  greatly 
except  at  the  20%  noise  level.  The  relative  error  for  recovery  of  r  (0-16%)  is  about  ten  times  greater 
than  the  relative  error  associated  with  the  recovery  of  es  (0-0.16%).  (This  agrees  with  relative  sensitivity 
of  the  model  to  the  parameters  seen  in  [10].) 
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Figure  17:  The  percent  relative  recovered  parameter  error  (100|r(true)  -r(recovered))/r(true)|)  as  a 
function  of  the  percent  noise  for  the  POD  (o,  50  modes)  and  FEM  (x)  formulations.  Data  points  are 
connected  by  straight  lines. 


Figure  18  plots  the  ratio  of  the  FEM  simulation  times  and  the  POD  simulation  times  as  a  function 
of  noise  (0-20%).  The  POD  formulation  is  4-10  times  faster  than  the  FE  method.  The  time  savings 
generally  decrease  as  the  noise  level  rises. 
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Figure  18:  The  ratio  of  solution  times  (Tfem/^pod)  as  a  function  of  the  percent  noise  for  the  two 
parameter  inverse  problem.  The  data  points  are  connected  by  straight  lines. 


4  REMARKS  AND  CONCLUSIONS 


The  saving  in  solution  time  using  the  reduced  order  POD  formulation  instead  of  the  FEM  formulation 
is  most  evident  as  the  number  of  nodes  increases.  The  solution  times  scales  as  N2,  where  N  is  the 
number  of  FEM  nodes  or  POD  modes.  However,  the  number  of  modes  necessary  to  capture  the  physical 
behavior  is  relatively  independent  of  the  number  of  nodes  used  to  generate  the  data. 

The  POD  method  is  significantly  faster  than  the  FEM  method  for  inverse  problems.  The  POD  method 
(50  modes)  is  4-10  times  faster  than  the  FEM  (149  nodes)  in  the  two  parameter  (es  and  r)  inverse 
problem  and  in  the  single  parameter  (es  or  r)  inverse  problem.  Based  on  our  results  from  the  forward 
simulation  problems  (where  savings  using  POD  over  FEM  dramatically  increased  as  N  increased,  e.g., 
50  fold  speed  up  for  N  =  450),  we  expect  the  speed  up  in  inverse  problems  to  dramatically  increase  in 
more  complex  problems  (e.g.,  2-D  and  3-D  problems)  where  an  increasing  number  of  basis  elements  are 
required. 

One  of  the  difficulties  with  using  the  POD  method  for  inverse  problems  is  in  the  construction  of  the 
POD  modes.  Memory  and  time  limitations  become  significant  as  the  number  of  parameters  is  increased 
in  a  time-dependent  problem  such  as  reported  here.  If,  in  a  IV-noded  simulation,  a  time  series  of  ( M ) 
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snapshots  is  generated  at  three  representative  values  spanning  the  desired  range  of  each  parameter,  there 
are  a  total  of  2>PMN  snapshots.  On  the  PC  used  for  these  calculations  there  is  a  memory  limitation  of 
approximately  6000  snapshots  (N  =  150)  in  the  algorithm  used  to  construct  the  POD  modes.  As  an 
example  of  the  time-limitiation,  construction  of  POD  modes  for  the  case  of  N  =  150  and  5010  snapshots 
required  18  hours  and  20  minutes  to  complete.  Current  efforts  and  ideas  to  alleviate  these  aspects  of 
difficulties  are  being  pursued  by  others  as  well  as  our  group  at  NCSU. 

The  findings  in  this  paper  offer  great  potential  and  encouragement  for  our  ongoing  efforts  in  higher 
dimensional  problems  where  a  number  of  significant  computational  difficulties  might  be  alleviated  by 
use  of  a  reduced  order  methodology  for  transient  electromagnetic  inverse  problems. 
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